Local equation of state and velocity distributions of a driven granular gas 
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We present event-driven simulations of a granular gas of inelastic hard disks with incomplete normal restitu- 
tion in two dimensions between vibrating walls (without gravity). We measure hydrodynamic quantities such as 
the stress tensor, density and temperature profiles, as well as velocity distributions. Relating the local pressure 
to the local temperature and local density, we construct a local constitutive equation. For strong inelasticities the 
local constitutive relation depends on global system parameters, like the volume fraction and the aspect ratio. 
For moderate inelasticities the constitutive relation is approximately independent of the system parameters and 
can hence be regarded as a local equation of state, even though the system is highly inhomogeneous with het- 
erogeneous temperature and density profiles arising as a consequence of the energy injection. Concerning the 
local velocity distributions we find that they do not scale with the square root of the local granular temperature. 
Moreover the high- velocity tails are different for the distribution of the x- and the y-component of the velocity, 
and even depend on the position in the sample, the global volume fraction, and the coefficient of restitution. 

PACS numbers: 45.70.-n, 51.30.+i, 51.10.+y 



I. INTRODUCTION 

The physics of vibro-fluidized granular materials is far from 
being fully understood. In particular, the applicability of hy- 
drodynamics 1 1, 2, 3, 4, 5] is still an object of debate |6, 7, 8]. 
Despite the similarities to the hydrodynamics of elastic hard- 
sphere systems, concerning e.g. the appearance of instabili- 
ties [l3 [ni Il2il . a main difference to ordinary fluids is 
the fact that continuous energy injection is vital to maintain 
a stationary state. Otherwise a gas of inelastic spheres would 
collapse, even in the absence of gravity. The non-trivial na- 
ture of this stationary state has been elucidated in experiments 
of vibrating grains. Among the most strikin g fe atures found 
are non-Gaussian vel ocity distributions lll3lll4l [Tsl fl6l| and 
cluster formation 11171 llSl Il9ll . Furthermore, several attempts 
have been made to extract either an equation of state or more 
gen erally a scaling relation for the thermodynamic variables 

The specification of the driving mechanism is a crucial in- 
gredient for any model describing vibrated granular fluids. A 
rather simple approach, though less appealing from an exper- 
imental point of view, utilizes stochastic bulk heating by un- 
corrected random forces, which act on every particle at every 
instant of time |23, S . In random restitution 

coefficients were considered with a probability distribution al- 
lowing for values both smaller and greater than one. However, 
this yields non-universal properties depending on the specific 
form of this distribution. The multiplicative bulk driving de- 
fined in terms of stochastic collision rules in JztIi is similar 
in spirit. As compared to bulk driving, a driving mechanism 
which acts only at a boundary of the system is much closer to 
experiments. In 101 the energy influx at the boundary has been 
modeled by a heuristically motivated ansatz for the heat cur- 
rent, while 1.28. 291 assume heating at the boundary through 
thermal walls. In the present paper the driving mechanism 
consists in incrementing a particle's velocity by a constant 
amount after every collision with a wall. This way of driving 
was used before in e.g. ifjcjlsHl . and with slight modifications 

in mill 13. 



One theoretical approach to study granular systems at fluid 
densities rests on kinetic theory, see e.g. d IH Ij? ] and 
references therein. Particular applications to driven granular 
gases can be found in |23, 2^ |27l |2^ El. Much of this 
work has been based on the Boltzmann equation, modified 
for inelastic collisions. One method to solve this nonlinear 
equation is based on the local equilibrium distribution 
which is only known for the elastic case. For systems with 
strongly inelastic collisions, the stationary state is unknown 
so that a systematic discussion of transport properties within 
kinetic theory is severely hampered. 

Hydrodynamic studies [HI Isoll have been moti- 
vated partly by the search for an understanding of temper- 
ature and density profiles |i2J, but also by experiments on 
hydrodynamic-like instabilities j^. In the hydrodynamic ap- 
proach the question arises, how to relate the pressure to the 
density and temperature. Several equations of state for the 
global quantities have been proposed, either interpolating be- 
tween the high and low density limit |2] or invoking the Boltz- 
mann equation lIsTIl for inelastically colliding particles. It is 
not clear, whether an equation of state also holds for the lo- 
cal hydrodynamic fields in a strongly driven, non-equilibrium 
system, where the fields are strongly inhomogeneous. 

Most simulations of vibrated granular gases have been 
based either on event-driven molecular dynamics L3Li40>.i4ll 
or the direct-simulation Monte Carlo method ll28ll30ll42ll43ll . 
If the system is driven through the boundaries, inhomoge- 
neous density and temperature profiles are measured. For low 
densities the computed temperature profiles agree well with 
hydrodynamic theory fl, "Su], whereas for moderate or high 
densities the profiles are not well understood with the excep- 
tion of almost elastically colliding particles. The full stress 
tensor, including potential contributions, has only been com- 
puted for freely cooling systems |44|. In simulations of driven 
granular gases |28] the collisional part of the stress tensor has 
not been measured, and in |31] the stress tensor is not mea- 
sured directly but instead computed from a local equation of 
state. 

Inspired by experiment 1 16], a lot of emphasis has been put 
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on the tails of the velocity distribution functions, which were 
found to be overpopulated as compared to a Maxwellian. In 
fact, all intermediate types of decay between a Gaussian and 
an exponential were observed llsil l42l ^45*. 46]. In addition 
mixtures [31. 47, 481 and rough spheres |49] have been inves- 
tigated as well as hydrodynamic instabilities such as convec- 
tion and pattern formation fs?] . 

In this paper, we present results from event-driven simu- 
lations of inelastic spheres in two dimensions confined be- 
tween two vibrating walls without gravity. Our focus is on 
the stationary state, which is reached when dissipation by 
particle collision equals energy injection due to the vibrating 
walls. In the stationary state, density and temperature profiles 
are shown to be strongly inhomogeneous due to the driving 
walls — even in a range of parameters, where clustering is only 
a minor effect. This has led us to 

(a) derive a constitutive equation by relating the mea- 
sured hydrodynamic fields, granular temperature T{x), vol- 
ume fraction (/)(a;) and pressure p{x), at each point x and 

(b) check whether this constitutive equation is universal or 
depends on the global system parameters of the model, like 
the aspect ratio of the cell, the overall volume fraction or the 
coefficient of restitution of the disks. For moderately inelastic 
systems (a = 0.9) the constitutive equation is (almost) inde- 
pendent of the remaining global system parameters so that the 
constitutive equation can be interpreted as a local equation of 
state for the a driven granular gas in the stationary state, even 
though the latter is highly inhomogeneous with heterogeneous 
temperature and density profiles. In contrast, for strongly in- 
elastic systems (a — 0.5) the constitutive equation depends 
significantly on the global system parameters so that the con- 
cept of a local equation of state cannot be sustained in this 
case. 

We furthermore discuss the one-particle distribution func- 
tion in the stationary state and show that 

(c) the local distribution fx{x,Vx) of Vx, the velocity in 
the direction of driving, is not a function of the rescaled vari- 
able Vxj \/Tx(x) alone. Similarly, curves of /y(a;, Vy) cannot 
be mapped onto a master curve for different x, when plotted 
against Vyj yjTy{x). Here Ti{x) denotes the local granular 
temperature associated with the translational motion in the i- 
direction. We find deviations from scaling at small and large 
arguments. 

(d) The local velocity distributions jx and have high- 
velocity tails whose decay ranges from stretched exponential 
to almost Gaussian. The particular type of decay depends on 
the position in the sample, the overall particle density and the 
coefficient of restitution. Furthermore, the decay of for 
large velocities is generally different from that of ]y. 

Taken together, one has to acknowledge that the stationary 
state of a driven granular gas is by no means universal. It 
shows peculiar features that depend on the precise values of 
the system parameters in contrast to the corresponding elastic 
system. 

The outline of the paper is as follows. In Sec. |ll|we in- 
troduce the model, specify the driving and define the observ- 
ables. We briefly discuss balance of energy input through 
the walls and energy dissipation in binary collisions in Sec. 
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FIG. 1: Model of N disks, driven in the x-direction with periodic 
boundary conditions in the y-direction. 



Unl Subsequently we present data for the profiles of the den- 
sity, the temperature and the components of the stress tensor 
(Sec. lIVI . In Sec.|V]we relate the local density, pressure and 
temperature to derive "experimentally" an equation of state 
and check for its universality. Finally, in Sec. lVII we discuss 
velocity distribution functions and their scaling behavior and 
present conclusions in Sec. lVIII 



II. MODEL AND OBSERVABLES 

We investigate a driven granular gas in 2 dimensions con- 
sisting of N identical inelastic smooth hard disks of diameter 
a and mass m which are confined to a rectangular box with 
edges of length Lx and Ly. The gas is driven through the walls 
perpendicular to the x-direction, which vibrate in an idealized 
saw-tooth manner (see below), while periodic boundary con- 
ditions are imposed in the y-direction. Fig.^shows a typical 
snapshot. The gas evolves in time through ballistic center- 
of-mass motion, binary inelastic collisions and particle-wall 
collisions. 



A. Binary collisions 

The inelastic nature of inter-particle collisions is the most 
important characteristic of granular media. As is often done 
we assume that it can be taken into account by a constant co- 
efficient of normal restitution and briefly recall the collision 
rules for binary collisions, see e.g. 

The unit-vector from the center of sphere two to the center 
of sphere one is denoted by fi := (ri — r2)/|ri — where 
Ti is the position vector of the center of mass of particle i. The 
center-of-mass velocities before a collision are denoted by Vi 
and V2, and the relative velocity by Vi2 = Vi — V2- Post- 
collisional quantities are primed. The relative velocity after a 
collision is assumed to obey 



n ■ V 



12 



-a n ■ Vi2 



(1) 



where a G [0, 1] is the (constant) coefficient of normal restitu- 
tion. The value a = 1 describes elastic collisions with energy 
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conservation, while for a < 1 energy will not be conserved 
but decreased in each collision. 

The constitutive equation ([0 plus momentum conservation 
determines the post-coUisional velocities for the disks 



v[ = Vi 



Av, 



pp' 



v'2 = V2 — At) 



pp 



where 



(l + a) . 
At>pp = (n • V12) n. 



(2) 



(3) 



The model could easily be extended to rough spheres by also 
specifying the disks' moment of inertia as well as a coefficient 
of tangential restitution /3o and Coulomb friction /i |55, 5^, 
but here we restrict ourselves to smooth spheres so that rota- 
tional motion is decoupled from translational motion. 



B. Driving 

When a particle collides with a driving wall, energy is in- 
jected into the system. This can be modeled in different ways, 
for example by drawing a new velocity from a Maxwellian 
distribution of a given wall temperature f?] or by assuming 
that the wall has a coefficient of normal restitution that is 
greater than one. Both of these mechanisms have no close 
experimental equivalents, though. A more realistic model is 
to assume a vibrating wall moving either in a symmetric (e.g. 
sinusoidal) or in an asymmetric (e.g. saw-tooth) way. In ad- 
dition, this can be combined with a normal (and also tangen- 
tial) coefficient of restitution |33, 41, 57]. In this article we 
refrain from the latter and restrict ourselves to saw-tooth driv- 
ing of the walls in the limit of vanishing amplitudes A and 
diverging frequency v such that Av =: Udrivc/2 is a constant. 
This ensures that the driving walls are always located at the 
same positions and leads to the following simple expression 
for a particle's change of velocity due to a collision with the 
left/right wall 

t)' = v + At>pw where Av^^ ^ {-2vx ±Vdnvc)ex (4) 
and B r stands for the unit vector in x-direction. 



C. Observables from simulations 



Following e.g. |I58|, 15911 we performed event-driven simu- 
lations in two dimensions with periodic boundary conditions 
in the ?/-direction and two identical idealized vibrating elastic 
walls in the x-direction. We initialize the system by placing 
the disks on a triangular lattice with a Gaussian velocity dis- 
tribution. To let the correlations of the initial state relax, the 
system is evolved elastically with periodic boundary condi- 
tions in the y-direction and elastic non-vibrating walls in the 
x-direction for an average of 100 collisions per particle. Then 
we switch on the driving of the left and right wall and the dis- 
sipation for particle-particle-collisions. Before we start mea- 
suring the observables we let the system relax further until, at 



time to, it has reached a stationary state, as indicated by the to- 
tal kinetic energy, which fluctuates around a time-independent 
mean value. 

To measure hydrodynamic fields in the stationary state we 
subdivide our box into cells Vr of area |Vr.|, centered at 
position r — {x,y). We then count the number of parti- 
cles f{r,Vx,Vy,t)\Vr\d,v.j;dvy at time t in cell Vr with an 
x-component of the velocity between and v.^ + dv^ and 
j/-component of the velocity between Vy and Vy + dvy. Such 
a local observable fluctuates as a function of time. To elimi- 
nate these fluctuations we average over a (long) time interval 
of length T and compute 

2 rto+r 

fstat{r,Vj;,Vy) -.^ - dt f{r,Vx,Vy,t) . (5) 

Of particular interest is the particle density 



(6) 



p{r) := / dvx / dvy fstat{r,Vx,Vy) 



or the area fraction <j>{r) := p{r)TTa?' /A and the two compo- 
nents i ~ x,y of the granular temperature 



p{r) 



dvx / dvy fst&t{r,Vx,Vy) - Vi{r)f 



where 



y(r) 



p{r) 



(7) 



dvx / dvy fsti,tir,Vx,Vy)v (8) 



is the velocity field. The total granular temperature in 2 di- 
mensions is defined as T{r) := {Tx{r) + Ty{r))/2. 

The stress tensor at position r and time t has a kinetic con- 
tribution and one that is due to the interactions between the 
particles cr(r, t) (T^™{r, t) + (T'"*(r, t). The kinetic part 
is given by 

o'-l'r(^) / d^x / dvy fstat{r,v.j;,Vy) 



x[v,-Vdr)][v,-V,{r)]. (9) 

If the particles interact through (finite) forces the contribu- 
tion due to interactions is given by the correlation of the par- 
ticles' relative positions and the forces between them. For 
hard-core interactions there are no forces, and one has to con- 
sider the momentum transfer in a small time interval instead 
L59i, i601: Suppose there is a collision at time t of particle k 
with another particle, then the momentum change of particle 
k will contribute to the stress tensor an amount proportional 
to (t) Apj (t). Here (t) is the i-th component of the vector 
of length a/ 2 pointing from the center of disk k to its col- 
lision contact point, and Ap'j{t) is the j-th component of the 
momentum change of particle k during this collision. To com- 
pute the collisional part of the stress tensor we need the total 
change of momentum in the time interval [t— At, t] in the cell 
Vr so that we have to keep track of all collisions n occurring 
at times t„ G [t — At, t], for which at least one collision part- 
ner kn (i.e. its center of mass) is located in cell Vr at time tn 
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^^^"^ = ip^EE^(^»)^p'"(*")- (10) 

The particle number fc„ of each such collision can take on one 
or two values, depending on whether one or both collision 
partners are located in cell V^. Similar to the other hydrody- 
namic fields, in a second step cr'"'(r, t) is averaged over time 
in order to get the collisional part of the stress tensor in the 
stationary state 

<J^(r)^lJ^''^\a^(r,t). (11) 

The corresponding local pressure p{r) —tr (cr(r)) /2 is 
defined as usual as the negative trace of the stress tensor di- 
vided by the space dimension. 

Besides hydrodynamic fields we have also measured veloc- 
ity distributions in Sec. IVII They are readily obtained by in- 
tegrating out the relevant variables in the stationary-state dis- 
tribution function /stat- A different method, which is better 
suited to determine high-velocity tails, for example, will be 
presented in Sec. lVII below. 

Coarse-grained measurements in space and time of certain 
observables may depend on the coarse-graining resolution. 
For example, this was demonstrated for the stress tensor in 
shear-flow driven granular systems in i60ll6lll . Our measure- 
ments of observables in the stationary state, however, should 
not suffer from such effects for two reasons: First, we could 
not detect any significant non-zero local velocity field in the 
simulated systems and, second, because of the long-time av- 
erage needed to obtain stationary-state quantities. 

Moreover, in the simulations we have never found any 
significant dependence of the (long-time-averaged) hydrody- 
namic fields on y, the coordinate parallel to the driving walls. 
Yet, stripe states, which are homogeneous in y, but have an en- 
hanced density in the middle of the sample, are known 1 10..1 1.1 
to exhibit instabilities with respect to density fluctuations in y. 
A marginal stability analysis 1 11 ] of granular hydrodynamics 
(for a close to one) determines the conditions under which 
such phenomena occur. As far as a comparison can be made, 
our systems fall into the stable region. Hence, we choose the 
cells Vr as stripes along the y-direction, for which we write 
Vx, and compute the hydrodynamic fields with spatial resolu- 
tion in x-direction only. Thus, we also write p{x) instead of 
p{r) and change the notation accordingly for all other quanti- 
ties. In the simulations the stripes Vx were all chosen to have 
equal width Lx/201, the temporal resolution of our measure- 
ments was set to At = 1, and the long-time average involves 
typically 10^ - 10^ collision events. 

Another instability of stripe states which is related to oscil- 
lations of the central dense cluster in x-direction |62| will be 
shortly addressed at the end of Subsec. lTV Al 

III. DIMENSIONAL ANALYSIS AND ENERGY BALANCE 

The model system contains three independent length scales: 
the diameter a of a disk and the box sizes Lx and Ly. In 




FIG. 2: Global granular temperature T as a function of the parameter 
defined below <12t . Comparison between simulations and the 
simple energy-balance argument, Eq. <12t . as well as with the refined 
version, Eq. <13> . The inset shows the same graph but on on a non- 
logarithmic scale. 



addition, there is one independent velocity scale, the driving 
velocity Udiivo, and one independent mass scale, the mass of 
a disk TO. Together with the initial positions and velocities of 
the particles that exhausts all dimensional quantities entering 
the time evolution of the system. 

We would like to describe the system using dimensionless 
variables that do not depend on the initial conditions since 
we expect stationary states to be independent thereof. There- 
fore, we will measure all lengths in units of the particle di- 
ameter a, all times in units of a/i;diivc and all energies in 
units of TOw^j.jy^. Note that there are no other time and en- 
ergy scales. Thus, we introduce dimensionless variables: box 
sizes Lx — Lx/a and Ly — Ly/a, granular temperatures 
Tx = Tx/{mvl^-^J and fy = Ty/{mvj^.^J, and the stress 
tensor cr = a a^/ (tow^j.^^^,). In the stationary state all dimen- 
sionless variables like T and a-, are independent of the driving 
velocity and only depend on position r and the remaining 4 
independent dimensionless system parameters, which charac- 
terize the system completely: the number of disks N , the two 
edges of length of the system Lx and Ly in units of a, and the 
coefficient of normal restitution a, which is a dimensionless 
material constant. 

For simplicity of the notation we refrain from indicating di- 
mensionless quantities by a tilde from now on. This should 
not cause confusion, because quantities having a physical di- 
mension will not occur any more in the rest of the paper (ex- 
cept for the appendix). 

We are interested in the macroscopic limit, which is taken 
such that N oo and Ly oo with a fixed line density 
A := N/Ly of particles. Thus, in the macroscopic limit the 
number of system parameters is further reduced to Lx, A, and 
a. It is important to keep Lx finite, otherwise energy bal- 
ance would not work: energy input occurs only at boundaries, 
while energy dissipation is a bulk phenomenon. 

It is instructive to estimate the average or global granular 
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temperature T := /_£ /2 P^^^ ^(^) balancing the en- 
ergy input at the walls and the energy loss due to particle col- 
lisions in the bulk |33, 63]. As shown in the appendix, we 
find that the granular temperature in the stationary state is in- 
dependent of the initial data and given by 

(l + v/l + (^2)2^)' (12) 

where ^jj := \/2x\{l — a^) and x stands for the pair correla- 
tion at contact of the corresponding elastic system, which we 
estimate by the Henderson approximation ( IA.6I 1. In Fig.|2]we 
plot the global granular temperature T as a function of ijj^^ 
and compare it with simulations. For small values of -0^^ 
the simulation deviates significantly from this simple theory's 
prediction. For large values of tA"^ (dilute and/or quasi-elastic 
systems) the agreement is reasonable. In addition, we show 
the result 

of a more refined calculation, which uses the pseudo- 
Liouville-operator approach to kinetic theory and will be pre- 
sented elsewhere. Eq. il3\ yields a better agreement with the 
data for intermediate values of tp~^. 

It is worth noting that the same type of argument which led 
to J12t also predicts that there is no generic stationary state 
for systems with the multiplicative driving mechanism (lA.U . 
described by Wdrivo = and a coefficient of restitution bigger 
than one for particle-wall collisions. Instead, such systems ei- 
ther cool down or heat up according to Haff's law i64ll . see 
Eq. ( IA.9I 1 in the appendix. This was also confirmed by simu- 
lations (not presented). 

IV. HYDRODYNAMIC FIELDS IN SIMULATIONS 

In this section we discuss the results of our simulations for 
the hydrodynamic fields, as computed from Eqs. (|6} - 
Due to the absence of a local velocity field V, these equations 
simplify accordingly. We have performed simulations for a 
wide range of system parameters a, (po, L^, Ly, and present 
examples thereof below. We then go on in Sec. [V] to discuss 
the question of a local equation of state, relating p{x), p{x), 
snAT{x). 

A. Density and temperature profiles 

We first present data for the density and temperature to 
demonstrate that the system is strongly inhomogeneous even 
for collisions with a = 0.9 and moderate densities. The 
parameters have been chosen as — 20, Ly = 25, 
and N = 256 so that the global area fraction is (po := 
TTa'^N/{AV) — 0.4. In Fig. |3lwe show the local area frac- 
tion </'(x), the X- and y-component of the granular tempera- 
ture, Tx{x) and Ty{x), as well as the isotropic temperature 
T{x) = [Tx{x) + Ty{x)]/2. We note that all quantities are 




FIG. 3: Spatial profiles of the area fraction (j> and granular tem- 
peratures Tx, Ty, and T from simulations [a — 0.9, Lx = 20, 
A = 10.24, A'^ — 256, (po = 0.4, corresponding to a mean-free path 
£/Lx ^ 0.02]. 



symmetric in x as expected. Except for a small area of ap- 
proximately one disk diameter next to the walls (indicated by 
the vertical lines), the area fraction is a monotonic function 
for either half of the system with the maximal value in the 
middle of the system. The temperatures T^, Ty, and T are 
monotonic as well with the lowest temperatures in the middle 
of the sample. 

The reason for the increased density next to the wall is an 
effective attractive potential of the wall due to entropic effects: 
Once a disks gets closer to the wall than one disk diameter, it 
can only receive hits from within the box but no hits from the 
direction of the wall. Thus the particle is pushed closer to 
the wall I65l l66ll . This effect is partially compensated by the 
driving walls, which add momentum to any particle hitting the 
walls. 

To support this explanation we have also investigated sys- 
tems of half the size —0.5 < x/L^ < 0, half the number 
of particles and with an elastic wall at a; = 0. The resulting 
hydrodynamic fields (not shown) are almost identical to the 
ones in Fig.|3]except for a small region of about one diameter 
close to .T = where the aforementioned effect is particularly 
visible. 

In systems exhibiting stripe states, we have sometimes ob- 
served an oscillatory instability of a central dense cluster in 
the x-direction |62], in particular if a is low and is large. 
However, these oscillations occur on far shorter time scales 
compared to the time interval r, over which we average to ob- 
tain stationary-state quantities. Hence, these oscillations are 
completely averaged out in the presented data. 



B. Stress tensor and pressure 

In Figs.|3-|6]we show the components of the stress tensor 
Cxxix), (Jyy{x), and (Jxy{x) for two different sets of parame- 
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FIG. 4: Spatial profiles of the stress tensor components a^,^, cjyy, FIG. 5: Same as Fig.|3 but for a less dense system [a = 0.9, = 
cr^y, and the product pT for the system of Fig.gl For all x away 50, A = 9.6, N = 240, 4>o = 0.15, l/L^ ^ 0.03]. 
from the boundary layers, pT is hardly to distinguish from pT^ and 
pTy (not shown). 

n — f ' 1 ' f ' f ' 1 — n 

0.0020 4- 



ters. Fig.l^is typical for systems with quite elastic collisions 
(a = 0.9) and moderate densities ((po = 0.4), whereas Figs.|5l 
and|6lshow rather dilute systems (<j>o = 0.15) at two different 
inelasticities, a moderate one (a — 0.9, as in Fig.|4} and a 
strong one (a — 0.5). The off-diagonal component cTj-j, of the 
stress tensor is always vanishingly small. The xx-component 
axx is constant within the sample except for a boundary layer 
close to the driving walls, whereas Fig.l^displays a dip of ayy 
in the center of the sample, which is more pronounced for the 
more dilute system in Fig.|5land hardly visible in the strongly 
inelastic, dilute system of Fig.|6l Furthermore, ayy increases 
considerably over a broad range in Fig.|6] when moving from 
a driving wall towards the center of the system. Currently, the 
origin of both x-dependences of (Tyy is not clear to us. While 
the stationary-state condition V • cr = and homogeneity in 
y require axx and axy to be constant in x on general grounds, 
this is not the case for ayy. We have carefully checked that 
the x-dependence of ayy is not caused by a shear instability 
associated with a non-zero velocity field V{x) in the system. 
Thus, there are normal stresses present in the simulated sys- 
tems, and they depend on x. 

If the equation of state of the ideal gas held, the local pres- 
sure p{x) would be related to the local temperature and den- 
sity according to p{x) = p{x)T{x) — {4/TT)(f){x)T{x). For 
our system the kinetic part of the stress tensor is diagonal and 
each component is simply related to the corresponding tem- 
perature a^i'ix) = pix)Tx{x) and a^l^'ix) = p{x)Ty{x). 
Hence the difference between the measured stress tensor and 
the ideal gas behavior is due to collisions. From Figs.|3-|6l 
it can be seen that the stress tensor deviates strongly from the 
ideal gas behavior. Consequently, collisions contribute signif- 
icantly. Here the collisional part has been measured directly 
in the simulations and not only estimated by approximate the- 
ories as has been done elsewhere ll30t 13111. 




We have also estimated the global mean free path £ for 



FIG. 6: Same as Fig.|5] but for a = 0.5. 



our simulated systems according to Eq. (9) in f?], which ex- 
presses £ solely as a function of the global volume fraction 
(f)Q. For the denser system with (/)o — 0.4, shown in Figs.|3l 
and0] this gives I k, 0.5 (in units of the diameter of the 
disks, that is t/Lx ~ 0.025). For the thinner systems with 
00 = 0.15, shown in Figs. |5] and |6j one gets I k, 1.65, that is 
l/Lx ~ 0.033. Thus, in all cases the mean free path £ is much 
smaller than the scales governing the spatial variations of the 
hydrodynamic fields, which are of order Lx- 



C. Density scaling 

For low global area fractions, 0o < 0.01, we observe scal- 
ing of the relative local area fraction (j>{x)/(j>o when plotted 
versus x/Lx- This is shown in Fig.0 where we compare den- 
sity profiles of systems with the same degree of inelasticity 
a = 0.9 and the same line density A — 10.24, but different 
values of the box width 1280 < Lx < 20000, correspond- 
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FIG. 7: Master curve for rescaled local area fraction in different low- 
density systems with the same line density [a — 0.9, A = 10.24, 
£/Lx ~ 0.03]. The dotted line corresponds to the theoretical predic- 
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FIG. 8: (Color online) Semilogarithmic parametric plots of the func- 
tion G(0) from the local equation of state <14> for different systems 
[a = 0.9, Ly = 25, 0.01 < ^/-L^ < 0.05]. Also shown are the 
theoretical predictions from <15> . <16> . and Eq. (4) in l67ll . 



ing to global area fraction 0o between 6 • 10~^ and 4 • 10~* . 
The relative local area fraction (/'/(/'o is seen to be a function 
of xjLx only and does not depend separately on and <^q. 
If the global area fraction were increased beyond 0o w 0.01 
(not shown), then the data collapse would cease to hold. In 
addition, the height of the peak would be reduced and instead 
of the bell-shaped master curve in Fig.0one would get overall 
concave profiles, similar to the one in Fig.|3l 

The hydrodynamic approach to quasi elastic, driven granu- 
lar gases in Sec. II D in |2] predicts for low-density systems 
that the mean free path in units of depends only on the 
Hne density according to IjL^ = 1/(V8A). When computed 
for the situation of Fig. one gets the small numerical value 
£/La: Ri 0.03, indicating that £ is not a relevant length scale 
for the master curve of the rescaled spatial density profiles. 
The theory of [2] also makes a prediction for the relative local 
area fraction (j){x) / 0o in terms of the solution of a first-order 
differential equation. This differential equation includes one 
free parameter, which we have fitted in order to match the so- 
lution (dotted line in Fig. at a; = with the master curve 
from our simulations. The agreement is reasonable, showing 
that inelasticities of a = 0.9 are at the borderline of the scope 
of this otherwise powerful approach to quasi elastic, driven 
granular gases. 



V. LOCAL EQUATION OF STATE 

Having measured hydrodynamic fields in the steady state 
of the system, one is led to search for a relation among them. 
Such a constitutive equation is needed, for example, in all hy- 
drodynamic approaches to driven granular gases in order to 
obtain a closed set of equations O lllll30ll . In the last section 
it was shown that a driven granular gas is intrinsically inho- 
mogeneous. Therefore it is only natural to investigate how the 
local values of the granular temperature T{x), pressure p{x), 



and density p{x) — ^resp. area fraction = p{x)tt/A in our 
units — are related to each other. To do so we observe in Fig.|3] 
that T{x), p{x), and (j){x) are all symmetric in x. Moreover, 
and this is crucial, (f){x) is monotone in x for either sign of x 
(except for a boundary layer of approximately one diameter 
in width close to a driving wall, which we ignore). Therefore 
one can invert the function (j>{x) for positive x. Upon inserting 
X = into the local pressure and temperature, we 

arrive at the constitutive equation 

with some function G. Fig. |8] shows parametric plots on a 
semilogarithmic scale of the function G with the values of 
p{x) / [p{x)T (x)] plotted against those of (p{x) for all x (ex- 
cept for those in the boundary layer mentioned above). In 
order to utilize a broad range of 0-values, Fig.|8lcontains data 
from 1 1 different systems, all of which have the same coef- 
ficient of restitution a = 0.9 and the same height Ly = 25, 
but different widths and different global area fractions (pQ. 
For not too large values of <j) these data merge quite nicely, 
indicating that there is only a weak dependence of G on the 
global system parameters and 0o in the corresponding pa- 
rameter range. In this case the constitutive equation ( I14> can 
be interpreted as the local equation of state of the system. The 
horizontal line in Fig.|8]marks ideal-gas behavior, from which 
G deviates due to the collisional contribution to the pressure. 
These deviations increase significantly with increasing (f>. 
The dotted line in Fig.|8lcorresponds to the function 

G(0) = l + + (15) 

where x stands for the pair correlation function at contact of 
the associated elastic hard-sphere gas in thermal equilibrium 
if^ . Since x is not known exactly, we estimate it by the Hen- 
derson approximation ( IA.6t for numerical purposes. In the 



8 



context of granular gases, il5\ occurred originally in a global 
equation of state for the ho mog eneous situation of a freely 
cooling granular gas Ql] . In i68ll the pair correlation function 
has been studied in a homogeneously driven inelastic system 
with periodic boundary conditions. It was found to be nearly 
independent of the coefficient of restitution and well approx- 
imated by the Henderson approximation, Eq. ( IA.6t . In |31] 
this form of G was used in the local equation of state (I14> 
to get a theoretical prediction for p{x) from simulations of 
T{x) and (j>{x) in a driven granular gas. Fig. |8] reveals that 
this works generally quite well for up to rather high local area 
fractions (/)(a;) < 0.5. In even denser systems agreement still 
holds for the well fluidized parts. Deviations from il5i start 
to occur when entering the transition zones to the frozen-out 
stripe of particles in the center of these high-density systems. 
The dashed-dotted line in Fig.|8]corresponds to the interpola- 
tion formula |2l 




for < (f> < (fic, which connects the behavior of dilute (van 
der Waals) and dense (ordered) elastic hard-sphere systems. 
Here 0c := 7r/(2V3) ~ 0.91 denotes the area fraction for 
ordered closed packings in two dimensions. Eq. ( I16t was ap- 
plied to quasi-elastic granular gases in |2], and even for our 
simulations with a = 0.9 there is agreement with the local 
data in the low-density regions up to (f){x) < 0.4. In addition, 
the most dense regions of regularly ordered, frozen-out par- 
ticles in the center of the high-density systems are described 
correctly, too. Yet another interpolation formula for G, which 
is rather accurate for an elastic hard-sphere gas even in the 
vicinity of the freezing transition, was put forward in Eq. (4) 
in |67] and is depicted by the dashed line in Fig.|S] Thus, the 
crossover between fluidized and frozen-out behavior in the in- 
elastic, driven systems, originating from the transition zone 
at the border of the solidified stripe of particles in the system 
center, is very smooth in comparison to the coexistence re- 
gion of the freezing transition for an equilibrated elastic gas. 
Moreover, the location of the crossover depends on the global 
system parameters. Since both interpolation formulae, il6\ 
and the one from f67'|, were tailored for elastic systems in 
thermal equilibrium, it is surprising to find such a good agree- 
ment with the local data for driven inelastic hard-sphere gases 
with restitution a = 0.9. 

When lowering the coefficient of restitution a from 0.9 to 
0.8 (not shown), there are only little changes to the plot in 
Fig- 13 In particular, the spread of the data pertaining to dif- 
ferent systems increases, even for low densities. Still one may 
be inclined to interpret il4l as a local equation of state — at 
least approximately. This situation is intermediate to the pre- 
vious with a = 0.9 and the following for strongly inelastic 
systems with coefficient of restitution a = 0.5. Indeed, Fig.|5] 
reveals major discrepancies in the local pressure, which are 
due to the systems with global area fractions above (po « 0.2. 
The discrepancies occur even at positions in the sample where 
the local area fractions are below 0.2 and where all the more 
dilute systems, i.e. those with 0o < 0.2, agree reasonably with 
the proposal ( I15t . The concept of a local equation of state is 
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therefore not sustainable any more for such strongly inelastic 
systems, and (I14> merely plays the role of a local constitu- 
tive equation, which depends in addition on the global system 
parameters. 

VI. ABSENCE OF SCALING FOR VELOCITY 
DISTRIBUTIONS 

Finally, we examine the local velocity distributions of the 
driven granular gas in the stationary state. We distinguish be- 
tween the local distribution 

f:^{x,Vx) —j—r [ dVy fstatix,Vx,Vy) (17) 

at position x of the velocity component Vx in the direction of 
the driving and the distribution 

fy{x,Vy) := ^— I dvx fstat{x,Vx,Vy) (18) 

P[X) JR 

of the velocity component Vy perpendicular to the direction 
of the driving. By definition, these velocity distributions are 
normalized to unity. 

In order to determine and fy from the simulation we use 
two different methods. The first one extracts them directly 
according to their definitions illl and ( I18> from /stat, which 
is determined as described in Subsec. Ill CI One disadvantage 
of this method is that it yields a smeared-out velocity distri- 
bution, which is spatially averaged over the width of the strip 
Vx centered around x. This is of practically no importance in 
the middle of the simulation box, but strongly disturbing for 
resolving the subtleties which occur close to the driving walls 
and are presented in Subsec. IVl Al The second way of mea- 
suring the velocity distributions avoids this problem: To get fi 
(i — xoT y) we keep track of all those particles which pass the 
line parallel to the y-axis at position x within a very long time 
interval of length r and whose ith component of the velocity 



lies in a small interval of width Av around Vi. Let us enumer- 
ate these particles by rii and denote their velocity components 



by V: 
Av 



(n) 

and Vy . Then one has (in the hmits r — > oo and 



0) 



p{x) fi{x,Vi) 



rAvL, 



(19) 



because for i ~ x, resp. 
equal to 



= y, the right-hand side of il9\ is 



\Fa:{x,Vx,Vy)\ 



resp. 



dVr 



\F^ix,Vx,Vy)\ 



(20) 

Here, F^{x,v.j,,Vy) := v^fstatix,v^,Vy) denotes the x- 
component of the (differential) current density in the station- 
ary state at position x of particles with velocity components 
Vx and Vy . 

As compared to the first method of measuring velocity dis- 
tributions, this one is also statistically more effective for de- 
termining rare events, such as the high-velocity tails in Suh- 
sec. lVlBl However, as far as fx{x, Vx) is concerned for small 
Vx, the second method is inferior to the first, because ( I19> 
assigns a large weight to the relevant events and therefore am- 
plifies statistical fluctuations, too. 



A. Effects of the discontinuity at a driving wall 

Particle-number conservation at a driving wall requires the 
incoming particle flux at the wall to be equal to the outgoing 
flux. For the velocity distribution fx this implies the boundary 
condition t30ll 



fxiTLx/2,Vx) =Q{±Vx - Wdr 



It 



Vdv 



X fxiTLxf^, -Vx + fdrivc) , (21) 

which must hold for all Vx > 0. Here Q{z) 1 for z > 0, re- 
spectively 8(z) := for z < 0, denotes Heaviside's unit-step 
function and Udrivc = 1 in our units chosen. The boundary 
condition i21\ relates the distribution fx of velocities prior 
to a collision with the wall to the one after a collision with 
the wall. In contrast, the velocity distribution fy must obey 
the usual reflection symmetry in Vy everywhere in the system, 
that is 



fy{x,Vy) — fyix, Vy) 



(22) 



and all Vy > 0. 



for all |a;| < 

We measured the velocity distribution fx{x, Vx) at 25 dif- 
ferent positions in a moderately inelastic system with coeffi- 
cient of restitution a = 0.9. Fig. 1101 shows the rescaled ve- 
locity distribution /^(x, w^;/ ^/TMaO) := fxix,Vx)\/Txix), 
measured using method one. At the driving walls fx is seen 
to obey the boundary condition ( I2U . When moving from a 
driving wall towards the center of the system, the gap in fx 
gets gradually smeared out, and fx becomes more and more 
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FIG. 10: Rescaled velocity distributions f^ for different x. [System 
parameters: a = 0.9, = 20, A = 10.24 (N = 256, (l>o = 0.4)] 




FIG. 11: 

left wall 



Two rescaled velocity distributions fx from Fig. 1101 at the 
and in the middle of the system. 



symmetric. Clearly, even for this moderately inelastic system 
fx does not scale for different x. The two extreme cases for 
X = -~Lx/2 and a; = are shown together in Fig.f^ The 
data for fx at the wall was obtained with the second method 
for measuring velocity distributions. This result is i n ag ree- 
ment with Direct-Simulation-Monte-Carlo results in 13011 and 
Molecular-Dynamics simulations in Isill . 



B. High-velocity tails 

Our main result for the velocity distributions is the non- 
scaling and multiformity of their tails. In order to observe 
these phenomena, extensive simulations for capturing rare 
events are required and the data have to be analyzed on a loga- 
rithmic scale. In contrast, on a linear scale the rescaled veloc- 
ity distributions seem to collapse approximately, as was ob- 
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FIG. 12: (Color online) Combined plot of reseated velocity distribu- 
tions of the same system as in Figs. 1101 and ll II fx in the middle of 
the sample (green crosses), fy at a driving wall (red squares), and 
fy in the middle (blue triangles). Also shown is a centered Gaus- 
sian with unit variance (black solid line). The inset shows the central 
part of the main graph with the (1433/2964/1340) data points being 
smoothed (running average over (41/85/38) data points). 



served previously by e.g. l69il . As an example we show 
in Fig. the rescaled distributions fy{x, Vy / ^jTy{x)) := 
fy{x, Vy)yjTy{x) in the middle of the sample (x = 0) and 
at a driving wall {x = ii^^/S) for the same moderately in- 
elastic system (a = 0.9) as in Figs. [TOJand^J Clearly, the 
curves are reflection-symmetric in Vy. For comparison, we 
have also included the symmetric distribution fx in the mid- 
dle of the system. Deviations are visible for small velocities 
only, as the inset of Fig.[T2lshows. The corresponding Gaus- 
sian (soUd line) in Fig. suggests that the velocity distribu- 
tions are close to but not identical to a Maxwellian. The ap- 
proximate data collapse, observed on this level of accuracy in 
Fig- El continues to hold if the coefficient of restitution 
of the gas is varied in a not too large extent. In considerably 
more inelastic systems, such as for a = 0.5, this is not true 
any longer For example, the peak of measured in the mid- 
dle of the system would be considerably broader and flatter 
than the ones of fy in the center and at the wall (not shown). 

In contrast. Fig. El ^) shows the same data of Fig. ll2lon 
a semi-logarithmic scale and includes also velocities of much 
higher absolute values. From this figure it is evident that scal- 
ing does not hold in the high-velocity tails of the distribu- 
tions either. Similar observations were made before in e.g. 
llsTl I45I l46ll . The type of decay in the high-velocity tails is 
different for fx and fy, and also depends on the position in 
the sample, the coefficient of restitution a, and the global area 
fraction 00 . This is illustrated by examples of different sys- 
tems in Fig. El In the insets we show | ln/i|~^ on a double- 
logarithmic scale in order to determine the decay exponent (3 
defined by 



ln/,(0,^;,/Vr,(0)) 



\Vy/^T^)\ 



(23) 



for fy in the middle of the sample. The exponent is defined 
accordingly for fy at the wall and f^ in the middle of the sam- 
ple. For the moderately inelastic systems with a = 0.9 in the 
first row of Fig. El the asymptotics has been clearly reached. 
We note that (3 is different for the different distributions, and 
also depends on the global area fraction 0o- Upon lowering 
a (top to bottom in Fig. El and/or decreasing the global area 
fraction 0o (left to right), the tails of fy get more and more 
populated, that is, /3 decreases. In very simple terms, this may 
be understood from the fact that (i) the largest typical veloci- 
ties are always of the order of f drive — 1 (see also Fig.ll4> and 
(ii) that Wdrivc/ \/Ti{x) increases up to 20 with decreasing a 
and decreasing (/)o. Hence, a Maxwellian velocity distribution 
would not be able to supply enough probability to particles 
with velocities of the order of I'diivc. instead higher-populated 
tails are needed. This argument suggests different behavior in 
different velocity regions so that the distributions cannot be 
fitted to the functional form j23t over the entire range of ve- 
locities 1 45]. Indeed, such a behavior can be seen in Fig.llBId). 
e) and f). The final asymptotics could not always be deduced 
from the simulations, even though our data include velocities 
which are up to 40 times bigger than the appropriate granular 
velocities ^/Ti. This applies to fy in the middle of the sam- 
ple in part e), where we suspect that the final asymptotics has 
not been reached. An asymptotic analysis of fx in the middle 
of the sample is even more problematic due to particles that 
reach the system center from a driving wall without undergo- 
ing a collision. These particles give rise to the side peaks of fx 
in parts e) and f), which have prevented us from determining 
f3 in these cases. 

Fig.ll3lcontains smoothed data as in the inset of Fig. El and 
the scales of the horizontal axes are determined by the square 
root of the appropriate granular temperatures. For compar- 
ison. Fig. El shows the unsmoothed data corresponding to 
Fig. Elf)' plotted directly versus Vi (in units of Wdiivo)- The 
side peaks of fx in the middle of the sample are due to the 
above-mentioned particles which fly from a driving wall to 
the system center without undergoing a collision. For the vast 
majority of particles |i>j^| is less than 80% of the driving ve- 
locity. It is in the region of this highest velocity observed for 
fy in the middle of the sample, where fx changes abruptly its 
slope. 

The exponent /3 has also been determined experimentally 
in a strongly driven gas so that gravity effects are small 1 16]. 
A value (3 w 1.55 was measured for a gas with coefficient 
of restitution a w 0.93. It was found to be remarkably in- 
dependent of the global area fraction (/)o, which was varied 
from 0.05 to 0.3. We have also simulated the system of llT6ll 
with zero gravity (not shown) and reproduced their value for 
/3. Even though the experimental data cover a wide range of 
velocities, some of the interesting phenomena discussed in 
Fig. El cannot be observed in this range. For the same rea- 
son the significance of the particular value /3 « 1.55 should 
not be overestimated. 

Simulations of a driven granular gas in a circularly shaped 
box also yield stretched Gaussian tails 1I45I1 . In addition, evi- 
dence is given that (3 depends only on the coefficient of resti- 
tution and the average ratio of the number of particle-wall 
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FIG. 13: (Color online) Rescaled velocity distributions fx in the middle of the sample (green crosses), fy at a driving wall (red squares), and 
fy in the middle (blue triangles) as in Fig. ^| but here smoothed data are shown on a semi-logarithmic scale and velocities of much higher 
absolute values are included. The solid line is a centered Gaussian with unit variance. Different parts of the figure represent different systems 
with decreasing global area fraction (po (left to right) and decreasing coefficient of restitution a (top to bottom). The remaining system data 
are = 20, A = 10.24 (iV = 256) for the left column and = 50, A = 9.6 (N = 240) for the right column. The insets show j In fi\~^ 
on a double-logarithmic scale to determine the decay exponent f} from <23> . The dashed line is the best linear fit to fx in the middle of the 
sample, the solid line to fy in the middle, and the dotted line to fy at the wall. Part a) corresponds to the system of Fig. fT2l 
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FIG. 14: (Color online) Velocity distributions for the system in 
Fig. 1131 0, but with unsmoothed data and without rescaling of the 
horizontal axis. 



to what extent the additional parameters affect the local equa- 
tion of state. In this context it would be interesting to study 
the effects of gravity for better comparison with experiments. 
Another interesting quantity to measure is the heat flux which, 
like the stress tensor, is also expected to be affected by colli- 
sions. Finally, one might also investigate other driving mech- 
anisms, like vibrations with finite frequency and non-zero am- 
plitude or sinusoidal driving. Work along these fines is in 
progress. 
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colUsions and particle-particle colUsions. For homogeneously 
driven systems the tail of the velocity distribution was the- 
oretically predicted |23] to be governed by the decay expo- 
nent /9 = 3/2. Simulations of homogeneously driven systems 
fTOl lTlll observed the exponent /? = 3/2 only for unrealisti- 
cally low values of the coefficient of restitution. 



VII. SUMMARY AND OUTLOOK 

In this article we explored the steady-state properties of a 
granular gas driven by vibrating walls. We have measured the 
full stress tensor, including the collisional contribution. This 
allowed us to obtain a constitutive equation directly from the 
simulations. The constitutive equation relates local density, 
pressure and temperature. For small inelasticities it can be 
regarded as the local equation of state of the gas, because it 
is to a large extent independent of the global system parame- 
ters. For strongly inelastic systems this interpretation cannot 
be sustained, instead the constitutive relation depends on the 
global volume fraction and sample geometry. We have also 
measured local velocity distributions, whose high-velocity 
tails were found to depend on the position in the sample as 
well as on the coefficient of restitution and the global volume 
fraction. Moreover, the tails are different for the two velocity 
components parallel and perpendicular to the driving walls. 
To conclude, the stationary state of a driven granular gas is in 
general non-universal — in contrast to the corresponding elas- 
tic system. This is not unexpected because the driven granular 
gas is a non-equilibrium state. Furthermore driving the sys- 
tem by energy input through the walls is effective only if the 
distance between the driving walls is finite so that the sample 
geometry enters naturally. 

We plan to extend our studies in various directions. It is 
straightforward to also consider rotational motion of the disks, 
by generalizing the collision rules to include tangential resti- 
tution as well as Coulomb friction. Again, the question arises 



APPENDIX: A SIMPLE ENERGY-BALANCE ARGUMENT 

Here we derive the energy balance Eq. il2\ by using argu- 
ments as in issl l63ll . We do this in a slightly more general 
setting than needed for Eq. il2\ and allow in addition for in- 
elastic collisions with the wall, characterized by a coefficient 
of restitution a„. The appropriate generalization of the colli- 
sion rule includes both the driving velocity and the coeffi- 
cient of restitution with the wall, 

v' = v + Avp^ with Avpy, = [-{l + a^)vx±Vdrivc]ex ■ 

(A.l) 

The special case Udiivo = and aw > 1 provides an alter- 
native driving mechanism, which, however, does not give rise 
to a stationary state, as will be shown below. In order to be 
able to treat this limit and for a better readability of the for- 
mulae, we refrain from using dimensionless physical units in 
this appendix. 

The average energy gain Ai?pw due to a particle-wall colli- 
sion is estimated from Eq. iA.l\ by averaging the kinetic en- 
ergy before and after the collision with a Maxwellian velocity 
distribution with (global) temperature T. This gives 



4av 



T ^ T 

^drivc \ ^ / 

ZTrm m 



(A.2) 

The collision frequency of particles with the left (right) wall 
is estimated by 



/p 



N 



T 
27rm ' 



(A.3) 



where we have assumed the density to be spatially homoge- 
neous throughout the system. 

When two disks collide in the bulk, the average change in 
total energy is computed similarly to ( IA.2> from (|2j and Q 
by averaging pre- and post-colhsional kinetic energy, which 
yields 

l-a2 



T . 



(A.4) 
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Finally, the number of particle-particle collisions per time is 
given approximately by Enskog's collision frequency I23il 



/pp = >^ax\l — : 



(A.5) 



where x is the pair correlation function at contact of a corre- 
sponding elastic gas in thermal equilibrium. Since x is not 
exactly known, we resort to the widely used Henderson ap- 
proximation I72I1 



X 



1 - 7(/)/16 



(1 



(A.6) 



for numerical purposes. It may be viewed as a heuristic ra- 
tional approximation to the virial expansion of x ™d is the 
two-dimensional equivalent to the Carnahan-Starling approx- 
imation 1661 for a three-dimensional hard-sphere gas. Ad- 
ditional higher-order terms to the Henderson approximation, 
which are proportional to (p^ / (1 — (p)^, have turned out to be 
irrelevant for our purposes and will therefore not be taken into 
account. 

Summing over energy loss in the bulk and energy gain at the 
right and left wall, we obtain for the total change in granular 
temperature 



dT 
It 



2 J'pw-^pw 



+ /pp-E'pp 



N 



m 



T 
2iTm 



+ 4q;^ 



T 
27rm 



2m 



(A.7) 



Here the dimensionless parameter 

2 



ipee ■=ip ("w - 1) 



(A.8) 



is given in terms of ip := \/2\ax{l — a^), which was already 
introduced below Eq. il2\ . (We note that Xa is the dimen- 
sionless line density employed there.) We briefly discuss two 
special cases: 

a) For Wdiivo = and a„ > 1 no stationary state is reached 
in general. Both, energy gain and loss increase like r3/2 

re- 
sulting in Haff 's law 

dt ~ 2L, V 2m ' ^ 

which has been discussed extensively in the different context 
of freely cooling granular gases. Here the temperature contin- 
ues to decrease or increase depending on whether dissipation 
or driving wins. 

b) For ctw = 1 and Wdrive > the granular temperature 
adjusts to the driving so that the stationary state with dT j dt = 
is characterized by the quadratic equation 



(A. 10) 




The solution of (|^^0} for the dimensionless global tempera- 
ture T :— T/ (mw^j.;^^) is given in Eq. ( fT2t . 
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